function [AB, AL, AC, AR, AT, S, dx, dy, dv, dc, siga, sct] = ...
    twoDcoefMG( dat, abdat, numg, numm, xcm, xfm, ycm, yfm, mt, src, it, BC, ...
    IBSL, IBSR, IBSB, IBST)
% function A = twoDcoef( dat, numg, numm, xcm, xfm, ycm, yfm, mt, src, it, BC,
%    IBSL, IBSR, IBSB, IBST)
%   This function produces the coefficient matrix and source for solution
%   of the two-dimensional diffusion equation.
tt = 0; ttt=0; ttv=0.000; tttv=0.000;
% expand the bc's
BCL = BC(1);
BCR = BC(2);
BCB = BC(3);
BCT = BC(4);

%---------------------------------------------------------------------temp
% determine legendre polynomial for incoming current
N = sum(xfm);    % max degree


% define data components
dc = zeros(numm,numg); sa=dc; ab = dc;
sc = zeros(numm,numg,numg);
for m = 1:numm
    for g = 1:numg %g+(m-1)*g
        dc(m,g) = dat(numg*(m-1)+g,1); % diffusion coefficient
        sa(m,g) = dat(numg*(m-1)+g,2); % absorption cross-section
        ab(m,g) = abdat(numg*(m-1)+g,1); % absorption cross-section
        ns(m,g) = dat(numg*(m-1)+g,3); % fission cross-section
        for gg = 1:numg
            sc(m,g,:) = dat(g+(m-1)*g,4:end);
        end
    end
end

% number of fine meshes in x and y coordinates
N = sum(xfm);
M = sum(yfm);
% number of coarse meshes in x and y coordinates
CN = length(xfm);
CM = length(yfm);
% compute dx and dy vectors
dx = zeros(N,1);
dy = zeros(M,1);
dv = zeros( (N+1)*(M+1), 1 );
j = 0;
for i = 1:CN
    dx( (j+1):(j+xfm(i)) ) = (xcm(i+1)-xcm(i))/xfm(i);
    j = sum(xfm(1:i));
end
j = 0;
for i = 1:CM
    dy( (j+1):(j+yfm(i)) ) = (ycm(i+1)-ycm(i))/yfm(i);
    j = sum(yfm(1:i));
end

lenk = (N+1)*(M+1);
AC = zeros( lenk, numg );
AL = zeros( lenk-1, numg );
AR = AL;
AB = zeros( lenk-N-1, numg );
AT = AB;
S = zeros( (N+1)*(M+1), numg ); siga = S;
sct = zeros( (N+1)*(M+1), numg, numg );

for g = 1:numg % group loop
    
    % interior coefficients
    for i = 2:N
        for j = 2:M
            for ci = 1:CN
                if i-1 <= sum(xfm(1:ci)), mtx = ci; break; end
            end
            for ci = 1:CM
                if j-1 <= sum(yfm(1:ci)), mty = ci; break; end
            end
            for ci = 1:CN
                if i  <= sum(xfm(1:ci)), mtxo = ci; break; end
            end
            for ci = 1:CM
                if j  <= sum(yfm(1:ci)), mtyu = ci; break; end
            end
            mtt  = mt(mtx,mty);   mtu  = mt(mtx,mtyu);
            mto  = mt(mtxo,mty);  mtb  = mt(mtxo,mtyu);
            
            k = i+(j-1)*( N+1 ); % 10 being the number of x ints
            
            alft = -0.5*( dc(mtt,g)*dy(j-1) + dc(mtu,g)*dy(j) ) / dx(i-1);
            argt = -0.5*( dc(mto,g)*dy(j-1) + dc(mtb,g)*dy(j) ) / dx(i);
            abot = -0.5*( dc(mtt,g)*dx(i-1) + dc(mto,g)*dx(i) ) / dy(j-1);
            atop = -0.5*( dc(mtu,g)*dx(i-1) + dc(mtb,g)*dx(i) ) / dy(j);
            a    = 0.25*( dx(i-1)*dy(j-1)*sa(mtt,g) + dx(i)*dy(j-1)*sa(mto,g) + ...
                dx(i-1)*dy(j)*sa(mtu,g) + dx(i)*dy(j)*sa(mtb,g) ) - (alft+argt+abot+atop);
            AL( k-1, g )   = alft;
            AR( k, g )     = argt;
            AB( k-N-1, g ) = abot;
            AT( k, g )     = atop;
            AC( k, g )     = a;  % ac ad  bc bd
            sct(k,g,:) =  0.25*( dx(i-1)*dy(j-1)*sc(mtt,g,:) + dx(i)*dy(j-1)*sc(mto,g,:) + ...
                dx(i-1)*dy(j)*sc(mtu,g,:) + dx(i)*dy(j)*sc(mtb,g,:) );
            if it == 0
                S( k, g )      =  0.25*( dx(i-1)*dy(j-1)*src(mtx,mty,g) + dx(i)*dy(j-1)*src(mtxo,mty,g) + ...
                    dx(i-1)*dy(j)*src(mtx,mtyu,g) + dx(i)*dy(j)*src(mtxo,mtyu,g) );
                dv( k, 1 ) = 0.25*( dx(i-1)*dy(j-1) + dx(i)*dy(j-1) + ...
                    dx(i-1)*dy(j) + dx(i)*dy(j) );
            else
                S( k, g )  = 0.25*( dx(i-1)*dy(j-1)*ns(mtt,g) + dx(i)*dy(j-1)*ns(mto,g) + ...
                    dx(i-1)*dy(j)*ns(mtu,g) + dx(i)*dy(j)*ns(mtb,g) );
                siga(k,g)  = 0.25*( dx(i-1)*dy(j-1)*ab(mtt,g) + dx(i)*dy(j-1)*ab(mto,g) + ...
                    dx(i-1)*dy(j)*ab(mtu,g) + dx(i)*dy(j)*ab(mtb,g) );
                if g == 1
                    dv( k, 1 ) = 0.25*( dx(i-1)*dy(j-1) + dx(i)*dy(j-1) + ...
                        dx(i-1)*dy(j) + dx(i)*dy(j) );
                end
            end
        end
    end
    
    % left and right edges
    for j = 2:M
        for ci = 1:CN
            if 2 <= sum(xfm(1:ci)), mtx = ci; break; end
        end
        for ci = 1:CM
            if j-1 <= sum(yfm(1:ci)), mty = ci; break; end
        end
        for ci = 1:CN
            if N-1  <= sum(xfm(1:ci)), mtxo = ci; break; end
        end
        for ci = 1:CM
            if j  <= sum(yfm(1:ci)), mtyu = ci; break; end
        end
        mtt  = mt(mtx,mty);   mtu  = mt(mtx,mtyu);                         %%%
        mto  = mt(mtxo,mty);  mtb  = mt(mtxo,mtyu);                        %%%
        i = 1;
        k = i+(j-1)*( N+1 );
        if BCL == 1
            AR( k, g ) = -1;
            AC( k, g ) = 1;
        else
            % average d in vert
            AR( k, g ) = -0.25 * ( dc(mt(1,mty),g) + dc(mt(1,mtyu),g) ) / dx(i);
            AC( k, g ) = 0.25 + 0.5 / dx(i) * ( dx(i)*dy(j-1)*dc(mt(1,mty),g) + dx(i+1)*dy(j-1)*dc(mt(mtx,mty),g) +...
                dx(i)*dy(j)*dc(mt(1,mtyu),g)  + dx(i+1)*dy(j)*dc(mt(1,mtyu),g) ) / ...
                ( dx(i)*dy(j-1)+dx(i+1)*dy(j-1)+dx(i)*dy(j)+dx(i+1)*dy(j) );
            if BCL == 2
                S( k, g ) = 0.5*(IBSL(mty,g)+IBSL(mtyu,g));
            end
        end
        
        S( k, g )  = tt*( 0 + dx(i)*dy(j-1)*ns(mto,g) + ...             %%%
            0 + dx(i)*dy(j)*ns(mtb,g) );
        siga(k,g)  = tt*( 0 + dx(i)*dy(j-1)*ab(mto,g) + ...             %%%
            0 + dx(i)*dy(j)*ab(mtb,g) );
        sct(k,g,:) = tt*(  dx(i)*dy(j-1)*sc(mto,g,:) + ...
            + dx(i)*dy(j)*sc(mtb,g,:) );
        if g == 1
            dv( k, 1 ) = ttv*( 0 + dx(i)*dy(j-1) + ...
                0 + dx(i)*dy(j) );
        end
        
        i = N+1;
        k = i+(j-1)*( N+1 );
        if BCR == 1
            AL( k-1, g ) = -1;
            AC( k, g )   = 1;
        else
            AL( k-1, g ) = -0.25 * ( dc(mt(end,mty),g) + dc(mt(end,mtyu),g) ) / dx(N);
            AC( k, g ) = 0.25 + 0.5 / dx(N) * ( dx(N)*dy(j-1)*dc(mt(end,mty),g) + dx(N-1)*dy(j-1)*dc(mt(mtxo,mty),g)  +...
                dx(N)*dy(j)*dc(mt(end,mtyu),g)  + dx(N-1)*dy(j)*dc(mt(mtxo,mtyu),g) ) / ...
                ( dx(N-1)*dy(j-1)+dx(N)*dy(j-1)+dx(N-1)*dy(j)+dx(N)*dy(j) );
            if BCR == 2
                S( k, g ) = 0.5*(IBSR(mty,g)+IBSR(mtyu,g));
            end
        end
        
        S( k, g )  = tt*( dx(i-1)*dy(j-1)*ns(mtt,g) + 0 + ...
            dx(i-1)*dy(j)*ns(mtu,g) + 0 );
        siga(k,g)  = tt*( dx(i-1)*dy(j-1)*ab(mtt,g) + 0 + ...
            dx(i-1)*dy(j)*ab(mtu,g) + 0 );
        sct(k,g,:) = tt*( dx(i-1)*dy(j-1)*sc(mtt,g,:) + ...
            dx(i-1)*dy(j)*sc(mtu,g,:) );
        if g == 1
            dv( k, 1 ) = ttv*( dx(i-1)*dy(j-1) + 0 + ...
                dx(i-1)*dy(j) + 0 );
        end
    end
    
    % top and bottom edges
    for i = 2:N
        for ci = 1:CN
            if i-1 <= sum(xfm(1:ci)), mtx = ci; break; end
        end
        for ci = 1:CM
            if 2 <= sum(yfm(1:ci)), mty = ci; break; end
        end
        for ci = 1:CN
            if i  <= sum(xfm(1:ci)), mtxo = ci; break; end
        end
        for ci = 1:CM
            if M-1  <= sum(yfm(1:ci)), mtyu = ci; break; end
        end
        mtt  = mt(mtx,mty);   mtu  = mt(mtx,mtyu);                         %%%
        mto  = mt(mtxo,mty);  mtb  = mt(mtxo,mtyu);                        %%%
        
        j = 1;
        k = i+(j-1)*( N+1 );
        if BCB == 1
            AT( k, g )     = -1;
            AC( k, g  )    = 1;
        else
            AT( k, g )     = -0.25 * ( dc(mt(mtx,1),g) + dc(mt(mtxo,1),g) ) / dy(1);
            AC( k, g )     = 0.25 + 0.5 / dy(1) * ...
                ( dx(i-1)*dy(j)  *dc(mt(mtx,1),g)    + dx(i)*dy(j)  *dc(mt(mtxo,1),g) + ...
                dx(i-1)*dy(j+1)*dc(mt(mtx,mty),g)  + dx(i)*dy(j+1)*dc(mt(mtxo,mty),g) ) / ...
                ( dx(i-1)*dy(j) + dx(i)*dy(j) + dx(i-1)*dy(j+1) + dx(i)*dy(j+1) );
            if BCB == 2
                S( k, g ) = 0.5*(IBSB(mtx,g)+IBSB(mtxo,g));
            end
        end
        
        S( k, g )  = tt*( 0 + 0 + ...
            dx(i-1)*dy(j)*ns(mtu,g) + dx(i)*dy(j)*ns(mtb,g) );
        siga(k,g)  = tt*( 0 + 0 + ...
            dx(i-1)*dy(j)*ab(mtu,g) + dx(i)*dy(j)*ab(mtb,g) );
        sct(k,g,:) = tt*(  dx(i-1)*dy(j)*sc(mtu,g,:) + dx(i)*dy(j)*sc(mtb,g,:) );
        if g == 1
            dv( k, 1 ) = ttv*( 0 + 0 + ...
                dx(i-1)*dy(j) + dx(i)*dy(j) );
        end
        
        j = M+1;
        k = i+(j-1)*( N+1 );
        if BCT == 1
            AB( k-N-1, g ) = -1;
            AC( k, g )     = 1;
        else
            AB( k-N-1, g ) = -0.25 * ( dc(mt(mtx,end),g) + dc(mt(mtxo,end),g) ) / dy(M);
            AC( k, g )     = 0.25 + 0.5 / dy(M) * ...
                ( dx(i-1)*dy(M-1)*dc(mt(mtx,mtyu),g) + dx(i)*dy(M-1)*dc(mt(mtxo,mtyu),g) + ...
                dx(i-1)*dy(M)  *dc(mt(mtx,end),g)  + dx(i)*dy(M)  *dc(mt(mtxo,end),g) ) / ...
                ( dx(i-1)*dy(M-1) + dx(i-1)*dy(M) + dx(i)*dy(M-1) + dx(i)*dy(M) );
            if BCT == 2
                S( k, g ) = 0.5*(IBST(mtx,g)+IBST(mtxo,g));
            end
        end
        
        S( k, g )  = tt*( dx(i-1)*dy(j-1)*ns(mtt,g) + dx(i)*dy(j-1)*ns(mto,g) + ...
            0 +0 );
        siga(k,g)  = tt*( dx(i-1)*dy(j-1)*ab(mtt,g) + dx(i)*dy(j-1)*ab(mto,g) + ...
            0 + 0 );
        sct(k,g,:) = tt*( dx(i-1)*dy(j-1)*sc(mtt,g,:) + dx(i)*dy(j-1)*sc(mto,g,:) );
        if g == 1
            dv( k, 1 ) = ttv*( dx(i-1)*dy(j-1) + dx(i)*dy(j-1) + ...
                0 + 0 );
        end
        
    end
    
    % finally, the corners
    for ci = 1:CN
        if i-1 <= sum(xfm(1:ci)), mtxL = ci; break; end  % picks out material for xfm(2)
    end
    for ci = 1:CM
        if 2 <= sum(yfm(1:ci)), mtyB = ci; break; end % picks out material for yfm(2)
    end
    for ci = 1:CN
        if i  <= sum(xfm(1:ci)), mtxR = ci; break; end % picks out material for xfm(end-1)
    end
    for ci = 1:CM
        if M-1  <= sum(yfm(1:ci)), mtyT = ci; break; end % picks out material for yfm(end-1)
    end
    %
    
    i=1; j=1; k = i+(j-1)*( N+1 );      % BOTTOM LEFT
    if (BCL==1&&BCB==1)
        AR( k, g )     = -1;
        AT( k, g )     = -1;
        AC( k, g )     = 2;
    elseif( BCL==1&&BCB==0 )
        AR( k, g )     = -1;
        AC( k, g )     = 1;
    elseif( BCL==1&&BCB==2 )
        AC( k, g )     = 0.25 + 0.5*( dc(mt(1,1),g)/dy(1) );
        AT( k, g )     = -0.25 / dy(1) * ( dc(mt(1,1),g) + dc(mt(mtyB,1),g) );
        S( k, g )      = IBSL(1,g);
    elseif(BCL==0&&BCB==1)
        AT( k, g )     = -1;
        AC( k, g )     = 1;
    elseif(BCL==2&&BCB==1)
        AC( k, g )     = 0.25 + 0.5*( dc(mt(1,1),g)/dx(1) );
        AR( k, g )     = -0.25 / dx(1) * ( dc(mt(1,1),g) + dc(mt(mtxL,1),g) );
        S( k, g )      = IBSL(1,g);
    else % default
        AR( k, g )     = -0.25 / dx(1) * ( dc(mt(1,1),g) + dc(mt(mtxL,1),g) );
        AT( k, g )     = -0.25 / dy(1) * ( dc(mt(1,1),g) + dc(mt(1,mtyB),g) );
        AC( k, g )     = 0.5 + 0.5*( dc(mt(1,1),g)/dx(1)+dc(mt(1,1),g)/dy(1) );
        if BCL == 2
            S( k, g )  = S( k, g )+IBSL(1,g);
        end
        if BCB == 2
            S( k, g )  = S( k, g )+IBSB(1,g);
        end
    end
    
    S( k, g )  = ttt*( dx(1)*dy(1)*ns(mt(1,1),g) );
    
    siga(k,g)  = ttt*( dx(1)*dy(1)*ab(mt(1,1),g) );
    sct(k,g,:) = ttt*( dx(1)*dy(1)*sc(mt(1,1),g,:) );
    if g == 1
        dv( k, 1 ) = tttv*( dx(1)*dy(1) );
    end
    
    i=1; j=M+1; k = i+(j-1)*( N+1 );    % TOP LEFT
    if (BCL==1&&BCT==1)
        AR( k, g )     = -1;
        AB( k-N-1, g ) = -1;
        AC( k, g )     = 2;
    elseif(BCL==1&&BCT==0)
        AR( k, g )     = -1;
        AC( k, g )     = 1;
    elseif( BCL==1&&BCT==2 )
        AC( k, g )     = 0.25 + 0.5*( dc(mt(1,end),g)/dy(1) );
        AB( k-N-1, g ) = -0.25 / dy(end) * ( dc(mt(1,end),g) + dc(mt(1,mtyT),g) ) ;
        S( k, g )      = IBSL(end,g);
    elseif (BCL==0&&BCT==1)
        AB( k-N-1, g ) = -1;
        AC( k, g )     = 1;
    elseif(BCL==2&&BCT==1)
        AC( k, g )     = 0.25 + 0.5*( dc(mt(1,end),g)/dx(1) );
        AR( k, g )     = -0.25 / dx(1) * ( dc(mt(1,end),g) + dc(mt(mtxL,end),g) );
        S( k, g )      = IBSL(1,g);
    else
        AR( k, g )     = -0.25 / dx(1) * ( dc(mt(1,end),g) + dc(mt(mtxL,end),g) );
        AB( k-N-1, g ) = -0.25 / dy(end) * ( dc(mt(1,end),g) + dc(mt(1,mtyT),g) );
        AC( k, g )     = 0.5 + 0.5*( dc(mt(1,end),g)/dx(1) + dc(mt(1,end),g)/dy(end) );
        if BCL == 2
            S( k, g )  = S( k, g )+IBSL(end,g);
        end
        if BCT == 2
            S( k, g )  = S( k, g )+IBST(1,g);
        end
    end
    
    S( k, g )  = ttt*( dx(1)*dy(end)*ns(mt(1,end),g) );
    sct(k,g,:) = ttt*( dx(1)*dy(end)*sc(mt(1,end),g,:) );
    siga(k,g)  = ttt*( dx(1)*dy(end)*ab(mt(1,end),g) );
    if g == 1
        dv( k, 1 ) = tttv*( dx(1)*dy(end) );
    end
    
    i=N+1; j=1;  k = i+(j-1)*( N+1 );   % BOTTOM RIGHT
    if (BCR==1&&BCB==1)
        AL( k-1, g )   = -1;
        AT( k, g )     = -1;
        AC( k, g )     = 2;
    elseif(BCR==1&&BCB==0)
        AL( k-1, g )   = -1;
        AC( k, g )     = 1;
    elseif(BCR==1&&BCB==2)
        AC( k, g )     = 0.25 + 0.5*( dc(mt(end,1),g)/dy(1) );
        AT( k, g )     = -0.25 / dy(1) * ( dc(mt(end,1),g) + dc(mt(end,mtyB),g) );
        S( k, g )      = IBSR(1,g);
    elseif(BCR==0&&BCB==1)
        AT( k, g )     = -1; %top
        AC( k, g )     = 1; %center
    elseif(BCR==2&&BCB==1)
        AC( k, g )     = 0.25 + 0.5*( dc(mt(end,1),g)/dx(end) );
        AL( k-1, g )   = -0.25 / dx(end) * ( dc(mt(end,1),g) + dc(mt(mtxR,1),g) ) ;
        S( k, g )      = IBSR(1,g);
    else
        AL( k-1, g )   = -0.25 / dx(end) * ( dc(mt(end,1),g) + dc(mt(mtxR,1),g) ) ;
        AT( k, g )     = -0.25 / dy(1) * ( dc(mt(end,1),g) + dc(mt(end,mtyB),g) );
        AC( k, g )     = 0.5 + 0.5*( dc(mt(end,1),g)/dx(end)+dc(mt(end,1),g)/dy(1) );
        if BCR == 2
            S( k, g )  = S( k, g )+IBSR(1,g);
        end
        if BCB == 2
            S( k, g )  = S( k, g )+IBSB(end,g);
        end
    end
    
    S( k, g )  = ttt*( dx(end)*dy(1)*ns(mt(end,1),g) );
    sct(k,g,:) = ttt*( dx(end)*dy(1)*sc(mt(end,1),g,:) );
    siga(k,g)  = ttt*( dx(end)*dy(1)*ab(mt(end,1),g) );
    if g == 1
        dv( k, 1 ) = tttv*( dx(end)*dy(1) );
    end
    
    i=N+1; j=M+1;  k = i+(j-1)*( N+1 );   % TOP RIGHT
    if (BCR==1&&BCT==1)
        AL( k-1, g )   = -1;
        AB( k-N-1, g ) = -1;
        AC( k, g )     = 2;
    elseif(BCR==1&&BCT==0)
        AL( k-1, g )   = -1;
        AC( k, g )     = 1;
    elseif(BCR==1&&BCT==2)
        AC( k, g )     = 0.25 + 0.5*( dc(mt(end,end),g)/dy(end) );
        AB( k-N-1, g ) = -0.25 / dy(end) * ( dc(mt(end,end),g) + dc(mt(end,mtyT),g) );
        S( k, g )      = IBSR(1,g);
    elseif (BCR==0&&BCT==1)
        AB( k-N-1, g ) = -1;
        AC( k, g )     = 1;
    elseif(BCR==2&&BCT==1)
        AC( k, g )     = 0.25 + 0.5*( dc(mt(end,end),g)/dx(end) );
        AL( k-1, g )   = -0.25 / dx(end) * ( dc(mt(end,end),g) + dc(mt(mtxL,end),g) );
        S( k, g ) = IBSR(1,g);
    else
        AL( k-1, g )   = -0.25 / dx(end) * ( dc(mt(end,end),g) + dc(mt(mtxR,end),g) );
        AB( k-N-1, g ) = -0.25 / dy(end) * ( dc(mt(end,end),g) + dc(mt(end,mtyT),g) );
        AC( k, g )     = 0.5 + 0.5*( dc(mt(end,end),g)/dx(end) + dc(mt(end,end),g)/dy(end) );
        if BCR == 2
            S( k, g )  = S( k, g )+IBSR(end,g);
        end
        if BCT == 2
            S( k, g )  = S( k, g )+IBST(end,g);
        end
    end
    
    S( k, g )  = ttt*( dx(end)*dy(end)*ns(mt(end,end),g) );
    sct(k,g,:) = ttt*( dx(end)*dy(end)*sc(mt(end,end),g,:) );
    siga(k,g)  = ttt*( dx(end)*dy(end)*ab(mt(end,end),g) );
    if g == 1
        dv( k, 1 ) = tttv*( dx(end)*dy(end) );
    end
    
    
end % group loop

